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solve the Radiative Transfer Equation in arbitrary one-, two- 
or three-dimensional optical environments. The optical envi- 
ronments are first divided into contiguous regions, or ele- 
ments, with Monte Carlo techniques then being employed to 
determine the optical response function of each type of 
element. The elements are combined, and the iterative 
relaxation techniques are used to determine simultaneously 
the radiance field on the boundary and throughout the 
interior of the modeled environment. This hybrid model is 
capable of providing estimates of the under- water light field 
needed to expedite inspection of ship hulls and port facili- 
ties. It is also capable of providing estimates of the subaerial 
light field for structured, absorbing or non-absorbing envi- 
ronments such as shadows of mountain ranges within and 
without absorption spectral bands such as water vapor or 

CO bands. 
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METHOD AND PROGRAM PRODUCT FOR 
DETERMINING A RADIANCE FIELD IN AN 
OPTICAL ENVIRONMENT 

RELATED APPLICATIONS 5 

This patent application claims priority to and all advan- 
tages of U.S. Provisional Patent Application 60/497,666, 
which was filed on Aug. 25, 2003. 

10 

STATEMENT OF GOVERNMENT INTEREST 

The work that led to this invention has been in supported 
in part by a grant from NASA, Grant Number NAS 5-31716. 
Thus, the United States Government may have certain rights 15 
to this invention. 

FIELD OF THE INVENTION 

The present invention generally relates to a method and 20 
program product for determining a radiance field in an 
optical environment. The present invention more specifi- 
cally relates to a method and program product that applies 
Monte Carlo techniques to determine the radiance field in 
the optical environment. 25 

BACKGROUND OF THE INVENTION 

In the many fields of physical science in which research- 
ers deal with aspects of radiative transfer, one important task 30 
is to determine the direction and intensity of radiant energy 
transport within a region of interest, given: 1 ) optical param- 
eters describing the scattering and absorption characteristics 
of the optical media within the region; 2 ) the specification of 
the direction and intensity of radiant energy incident on the 35 
boundary of the region; and 3) the specification of any 
radiant source distribution within the region. For example, 
the development of optical coherence tomography tech- 
niques in the realm of medical optics requires knowledge of 
the response of new optical instruments. 40 

Accurate calibration of blackbody sources used in pyrom- 
etry and spectroscopy requires determination of the emission 
characteristics of new experimental sources. Remote sensing 
of the earth’s atmosphere, terrain and marine environments 
depends upon instruments, the design and calibration of 45 
which are possible only after the specification of the optical 
conditions in which the instruments are to perform. These 
specifications are frequently determined by simulation of the 
radiant energy transport in the physical conditions to be 
encountered by the instruments. 50 

The process common to all three endeavors mentioned 
above, namely, the radiative transfer of energy through 
media capable of interacting with and altering the radiance 
field, is governed by the Radiative Transfer Equation or 
RTE. Unfortunately, the RTE supplies no analytic solution 55 
for problems characterized by realistic geometries, boundary 
conditions and optical parameters. As a result, the “solving” 
of the RTE for realistic optical environments consists cur- 
rently of specifying the parameters necessary to define the 
problem, and then numerically simulating the solution, such 60 
as the radiance field, at the desired location(s). 

Many methods for generating numerical solutions to the 
RTE in the optically coupled ocean-atmosphere environ- 
ment have been developed in support of studies of the 
earth’s atmosphere and oceans. Most of this environment 65 
can be represented by one-dimensional models that treat the 
atmosphere and oceans as stratified media that are horizon- 
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tally infinite in extent and in which the optical parameters 
vary only with depth or height. Indeed, in the common case 
of clear sky over open ocean waters, the one-dimensional 
approximation is very good; the one-dimensional models 
execute quickly and have negligible or acceptable error. In 
order to operate effectively, a researcher or other person 
using monitoring equipment should be capable of establish- 
ing parameters for the particular environment and the con- 
ditions at the time of taking readings. Because these param- 
eters are dynamic in nature, it is necessary to formulate a 
scenario for the conditions at the time of the instrumental 
measurements to determine the effectiveness of the results as 
well as the set up conditions for the instrumentation. This 
has been found to be possible in difficult media by using the 
Monte Carlo method which can give an accurate depiction 
of the magnitude of the measured entity for a given set of 
conditions. However, Monte Carlo calculations are also long 
and cumbersome in nature and require a great amount of 
data manipulation. 

SUMMARY OF THE INVENTION AND 
ADVANTAGES 

A method and program product for determining a radiance 
field in an optical environment are disclosed. The method 
includes generating a model of the optical environment and 
partitioning the modeled, optical environment into discreet, 
contiguous regions. The method further includes applying 
Monte Carlo techniques to determine an optical response 
function for each region for which the optical response 
function cannot be analytically defined. The program prod- 
uct includes codes for accomplishing the various steps 
associated with the method described above. 

The method and program product of the present invention 
utilize Monte Carlo techniques to solve a Radiative Transfer 
Equation (RTE) in arbitrary one, two, or three-dimensional 
environments. The present invention permits the rapid deter- 
mination of characteristics of an optical environment for a 
known and thus standard set of parameters. This enables a 
more efficient and finer processing of information to occur. 
From this information, it is possible to gain a perspective on 
the accuracy and actual ability of an optical instrument, such 
as a sensor, to perform an analysis for given conditions at 
any specific time. 

More specifically, the method and program product of the 
present invention approximate the solution of the RTE in 
any dimension. To effect this approximation, the method and 
program product generate a model of the optical environ- 
ment and partition the modeled, optical environment into 
discreet, contiguous regions and then solves the RTE in each 
element using Monte Carlo techniques. The output from 
each element then serves as the input to adjacent elements. 
Radiance from boundary and internal sources is diffused 
throughout the environment as the element outputs are 
readjusted iteratively. Once the iterative technique con- 
verges, the radiance is known on the boundary and through- 
out the entire modeled environment. The method of the 
present invention predicts the efficacy of instrumental mea- 
surements for a given set of conditions and provides a set of 
parameters for given situations which may be used to predict 
the outcome and error probabilities for those conditions. 

BRIEF DESCRIPTION OF THE SEVERAL 
VIEWS OF THE DRAWINGS 

Having thus described the invention in general terms, 
reference will now be made to the accompanying drawings, 
which are not necessarily drawn to scale, and wherein: 
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FIG. 1 depicts the transmission of irradiance into and out 
of a region of interest, according to an illustrative example 
of the present invention. 

FIG. 2 shows a 3-D cubic structure with each face 
subdivided into 16 grid cells, according to an illustrative 
example of the present invention. 

FIG. 3 is a 2-D plate, according to an illustrative example 
of the present invention. 

FIG. 4 shows the modeling of one dimensional environ- 
ments using three-dimensional cubes, according to one 
aspect of the present invention. 

FIG. 5 shows the modeling of an infinitely long bar of 
square cross-section, according to one aspect of the present 
invention. 

FIG. 6 shows a model geometry for solving a 1 -dimen- 
sional problem using long bars and strips, according to one 
aspect of the present invention. 

FIG. 7 shows a model geometry for solving a 2 -dimen- 
sional problem, according to one aspect of the present 
invention. 

FIG. 8 illustrates the downward irradiance (E rf ) vs. depth 
as calculated by the methods of the present invention and by 
Hydrolight, according to one illustrative example. 

FIG. 9 illustrates the upward Radiance (E w ) vs. depth as 
calculated by the methods of the present invention and by 
Hydrolight, according to one illustrative example. 

FIG. 10 shows the average cosines plotted vs. depth as 
calculated by the methods of the present invention and by 
Hydrolight, according to one illustrative example. 

FIG. 11 shows the total upward radiance just above the 
surface in the half-plane perpendicular to the principal plane 
as calculated by the methods of the present invention and by 
Hydrolight, according to one illustrative example. 

FIG. 12 shows an illustrative normalized Log 10 E rf at 532 
nm, Chi =0.3 mg/m 3 as calculated by the methods of the 
present invention. 

FIG. 13 shows another illustrative normalized Log 10 E w at 
532 nm, Chi =0.3 mg/m 3 as calculated by the methods of the 
present invention. 

FIG. 14 shows an illustrative normalized Log 10 E rf at 532 
nm, Chi =5.0 mg/M 3 as calculated by the methods of the 
present invention. 

FIG. 15 shows an another illustrative normalized Log 10 E w 
at 532 nm, Chl=5.0 mg/m 3 as calculated by the methods of 
the present invention. 

DETAILED DESCRIPTION OF THE 
INVENTION 

The present inventions now will be described more fully 
hereinafter with reference to the accompanying drawings, in 
which some, but not all embodiments of the invention are 
shown. Indeed, these inventions may be embodied in many 
different forms and should not be construed as limited to the 
embodiments set forth herein; rather, these embodiments are 
provided so that this disclosure will satisfy applicable legal 
requirements. Like numbers refer to like elements through- 
out. 

The present invention includes a method and program 
product for determining a radiance field in an optical envi- 
ronment. Generally, the method and program product gen- 
erate a model of the optical environment, partition the 
modeled, optical environment into discreet, contiguous 
regions, or elements, and apply Monte Carlo techniques to 
determine an optical response function for each region for 
which the optical response function cannot be analytically 
defined. The optical response function provides an output for 
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any given input. The step of applying Monte Carlo tech- 
niques to each region more specifically comprises solving 
the Radiative Transfer Equation (RTE) for each region to 
obtain an end solution indicative of the radiance field in the 
5 optical environment. 

Those skilled in the art readily recognize that there is not 
‘one’ RTE. Instead, there are many different forms of the 
RTE but, generally, the RTE is the equation which describes 
the propagation of radiation (e.g. light) through a scattering 
10 and absorbing medium. Radiation and light are used inter- 
changeably herein. However, it is to be understood that light 
is but one form of radiation and, in a general context, 
radiation refers to any type of electromagnetic radiation. 

15 For the purposes of the present invention, the terminology 
optical environment is directed at any environment, or 
region of space, that is desirable to model where particular 
attention is paid to optical properties of the environment 
including, but not limited to, scattering and absorption 
20 characteristics. Non-limiting examples include the optical 
environment underwater or subaerial. Further, it is under- 
stood by those skilled in the art that the optical response for 
some regions can be analytically defined. As such, applica- 
tion of the Monte Carlo techniques is not required toward 
25 these particular regions. In other words, the Monte Carlo 
techniques do not have to be applied to each region that has 
been partitioned. 

The method and program product further include com- 
bining the optical response function for each region. Further, 
30 the output of one region is utilized as an input to a contigu- 
ous region. The output of one region, therefore, depends on 
the input into the same region, as well as the optical response 
function. Preferably, iterative relaxation techniques are then 
applied to determine the radiance field on a boundary and 
35 throughout the modeled, optical environment. In this con- 
text, the outputs from each region are iteratively adjusted. 
Due to the application of the iterative relaxation techniques, 
the present invention is a hybrid method that includes two 
primary phases, (1) application of the Monte Carlo tech- 
40 niques followed by (2) application of the iterative relaxation 
techniques. The program product of the present invention 
includes specific codes, such as a modeling code, a parti- 
tioning code, an application code, a combination code, and 
a relaxation code for accomplishing the various steps asso- 
45 ciated with the method described above. 

The method of the present invention includes many 
different practical applications in response to obtaining the 
end solution. For instance, the method may calibrate data 
from an optical sensor in response to the end solution. As 
50 described additionally below, the calibration of the data may 
include adjusting gain of the optical sensor, preferably 
automatically, in response to the end solution. The method 
may validate data from an optical sensor in response to the 
end solution. Further, the method may determine a particular 
55 optical sensor to deploy in response to the end solution. 
Even further, the method may adjust an optical sensor in 
response to the end solution prior to deployment of the 
optical sensor. As alluded to above, the term optical sensor 
is used interchangeably with optical instrument. In any 
60 event, although not required, it is preferred that the optical 
sensor is located on at least one of a remotely operated 
vehicle (ROV), an autonomous underwater vehicle (AUV), 
a piloted aircraft, an autonomous aircraft, a semi-autono- 
mous aircraft, a piloted spacecraft, an autonomous space- 
65 craft, and a semi -autonomous spacecraft. Further, it is to be 
understood that the terminology spacecraft includes satel- 
lites. 
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In the field of instrumental analysis, many calibration 
methods are based on one-dimensional models. These mod- 
els fail, since in the case of field observations, the world is 
not one-dimensional in nature. In fact, the efficiency of the 
instrument is dependent of the environmental parameters 
encountered in the field. Some examples of this include the 
radiance received at a spaceborne sensor viewing coastal 
areas where relatively dark water adjoins bright beaches, or 
viewing maritime scenes in the presence of brightly reflect- 
ing clouds. One-dimensional models also fail when the 
geometry of the region to be modeled is not one-dimensional 
in nature, such as when modeling the light field adjacent to 
and beneath a ship. 

Monte Carlo techniques comprise the only known method 
of solving the RTE for those problems not amenable to the 
one-dimensional approximation. The radiance distribution 
in any region, including those with irregular geometries, 
parameter distributions, and boundary conditions, may be 
modeled by Monte Carlo methods. Furthermore, the solu- 
tion to the RTE in any region may be achieved to any desired 
degree of accuracy and may also be collated for a prescribed 
set of conditions. Thus Monte Carlo packages or modules 
for differing conditions may be formulated for those condi- 
tions and then may be used to predict the efficiency of 
instrumentation or the reliability of results for any set of 
conditions. In this way, an operational package or Monte 
Carlo code may be produced for the study conditions for a 
set of circumstances. 

In general, there are three serious disadvantages of Monte 
Carlo methods. First, for the most part, Monte Carlo models 
do not execute quickly. Monte Carlo techniques require the 
accumulation of statistics based on tracing the histories of 
vast numbers of photons from their source to their extinction 
or lack of relevance. The process involves repeated sampling 
of probability distributions to determine the type and out- 
come of events the photons undergo. In realistic three- 
dimensional regions the ray-tracing routines necessary to 
calculate the path of the photons, their intersections with 
objects and boundaries, and their reflections and refractions, 
can become computationally burdensome. Therefore, Monte 
Carlo codes may require hours, days, or even weeks to 
execute, depending on geometric complexity and the quan- 
tity of interest. 

Second, Monte Carlo codes are generally problem spe- 
cific in nature. Although some subroutines, such as those 
determining path length or scattering angle, may be reus- 
able, most “new” problems require extensive coding and 
testing. Once the new code compiles and executes correctly, 
it may be difficult to verily the model results, since, in 
theory, the model is simulating a scenario for which no other 
model exists and no physical measurements are available. 

Third, Monte Carlo models generally provide results for 
only one set of optical parameters at a time, and at once, or 
at most, a few, locations within a modeled region. To change 
optical parameters or to move the location of the model 
sensor system usually requires rerunning of the model. 

For these reasons, at least within the context of atmo- 
sphere/ocean remote sensing, Monte Carlo techniques are 
most often employed to solve the RTE in two- or three- 
dimensional situations for which the faster techniques are 
not applicable. Monte Carlo solutions for these cases may be 
used to determine the accuracy of one-dimensional models 
when they are applied to two- or three-dimensional sce- 
narios. Once the errors are quantified, corrections to one- 
dimensional models may be developed that allow the models 
to be applied to environments that are perturbed from 
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one-dimensional in specific ways, such as for sloped or 
patched sea floors underlying an otherwise one-dimensional 
water column. 

Methods and computer program products of the present 
5 invention approximate the solution of the RTE in any 
dimension. To effect this approximation, the methods and 
computer program products of the present invention parti- 
tion a model environment into contiguous regions, or ele- 
ments, and solve the RTE in each element using Monte 
to Carlo techniques. The output from each element then serves 
as the input to adjacent elements. Radiance from boundary 
and internal sources is diffused throughout the environment 
as the element outputs are readjusted iteratively. Once the 
iterative technique converges, the radiance is known on the 
15 boundary and throughout the entire modeled environment. 

The advantages of the instant method stem from the reuse 
of basic environmental elements, and from the ease with 
which the elements may be combined to form new geom- 
etries. Furthermore, once a particular environment is mod- 
20 eled, the time required to solve for a similar environment is 
substantially reduced. 

It will be appreciated that the methods of the present 
invention are described below with reference to block dia- 
grams. It will be understood that each method described 
25 below can be implemented by computer program instruc- 
tions. These computer program instructions may be loaded 
onto one or more general purpose computers, special pur- 
pose computers, or other programmable data processing 
apparatus to produce machines, such that the instructions 
30 which execute on the computers or other programmable data 
processing apparatus create means for implementing the 
functions specified in the flowchart block or blocks. Such 
computer program instructions may also be stored in a 
computer-readable memory that can direct a computer or 
35 other programmable data processing apparatus to function in 
a particular manner, such that the instructions stored in the 
computer-readable memory produce an article of manufac- 
ture including instruction means that implement the function 
specified in the flowchart block or blocks. Therefore, the 
40 methods described herein may be implemented by a mod- 
eling tool comprising a processor, operating system, 
memory, input/output interface, one or more databases and 
a bus. 

Further, the methods described herein may be embodied 
45 as a data processing computer-readable program code means 
embodied in the storage medium. Any suitable computer- 
readable storage medium may be utilized including hard 
disks, CD-ROMs, DVDs, optical storage devices, or mag- 
netic storage devices. Accordingly, the present invention 
50 may take the form of an entirely hardware embodiment, an 
entirely software embodiment or an embodiment combining 
software and hardware aspects, such as firmware. According 
to a preferred embodiment, die methods described herein are 
implemented by a stand-alone software application operat- 
55 ing on a Windows® operating system. However, the meth- 
ods may be implemented using alternative operating systems 
and databases as are known to those of skill in the art. 

Referring now to FIG. 1, methods of the present invention 
consider a point, 10, on the surface of an arbitrary region. 
60 This point 10, is illuminated by a beam of radiant power, 11, 
at wavelength 2 and from direction 12. E. is that portion of 
the iiradiance at point 10 that arrives from direction 12 and 
is expressed in units W m -2 sr -3 . As a result of this illumi- 
nation, power may be emitted from the region at 0 point 15, 
65 along direction 17 toward detector 16. 

The linearity of the RTE implies that the amount of energy 
leaving the region at point 15 and in direction 17 resulting 
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from the energy at point 10 from direction 12 is proportional 
to the amount of energy at point 10 from direction 12. The 
point-by -point constant of proportionality is denoted as the 
“response function” and is defined as: 

£ o (p o ,e^£/p z -,e z 0)=^(p o ,G o ;Pz,e^z(Pz,Gz)- 

If R is known for all points on the boundary of the region, 
and for all directions, then the power emitted by the region 
is completely determined by the incident power distribution. 
The total power emitted from differential of the surface area 10 
centered on point 15 and propagating within the differential 
cone of directions centered on direction 16 is given by: 

£ 0 (p 0 ,eJ^(pJ^Q(eJ=JJi?(p 0 ,e 0 ;p z .,e z .)£: z .(p z .,e z .)^ 

dQ 

The integrals in the above equation are expressed over the 1 5 
entire boundary of the region and all directions, dQ, which 
are incident on the region at d A. 

The relationships given in the above equations are useful 
for Monte Carlo techniques and have two properties. First, 

R is a relationship between points on the surface of the 20 
region and is not explicitly dependent on the optical media 
within. R depends on the material within the region. But 
when R can be measured, deduced, or defined for the region, 
when the power distribution incident on the surface of the 
region is known, the emitted power distribution can be 25 
calculated without further reference to the optical properties 
of the region. Secondly, there are regions for which R is 
known or may be defined exactly. These include regions 
within which no scattering takes place and regions with 
boundaries that have known reflectance characteristics. 

Constructing A Model 

The first step in creating a simulation of the light field in 
a given optical environment is to partition the environment 
into discreet, contiguous elements in a similar manner as is 35 
employed in a finite element analysis. Two element shapes 
may be used to illustrate the methods of the present inven- 
tion: the three-dimensional cube as shown in FIG. 2 and the 
two-dimensional square plate depicted in FIG. 3 . It will be 
appreciated that their shapes are for illustrative purposes 40 
only, and are non-limiting examples. 

In the illustrative example shown in FIG. 2, the 3-D cube 
has a directional resolution of 5°x5° with an input at face 2, 
cell 4, theta bin 24 and phi bin 68 . Generally, a cube may be 
used to represent a three-dimensional region of the model 45 
environment containing a well-defined distribution of opti- 
cally reactive constituents such as seawater, particulates, etc. 
Once the RF is calculated for this particular distribution of 
optical agents, that RF may be used repeatedly for other 
cubic regions of the model environment that contain the 50 
same constituent distributions. 

The illustrative plate shown in FIG. 3 represents a Fresnel 
interface with the reflected (E^) and transmitted (E r ) irra- 
diance results from an input E.. The plate may be used to 
represent, e.g., boundaries, interfaces, and sources. In this 55 
way, one or more plates may be used to represent the sea 
floor, the bottom and sides of barges, and/or an air/sea 
interface. Plates may also be used to represent the upper 
boundary of a modeled environment, such as a source with 
a radiance distribution representing the solar beam and 60 
skylight. 

The spatial resolution of a model depends on the resolu- 
tion of the special grid imposed on the surface of the 
elements as depicted in FIGS. 2 and 3. Each output value of 
the element is an average over the area of the grid with 65 
which it is associated. Therefore, it is expected that the grid 
resolution should be small compared to the scale of spatial 
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change of the expected modeled light field. For instance, a 
smaller grid size may be required to model the area near the 
air/sea interface than is required to model a region in deep, 
homogeneous water. In addition the size of the model (e.g., 
a chosen cube structure such as that shown in FIG. 2) may 
be selected to replicate sizes of the physical shapes found in 
the proposed model environment. As in FEM where a 
problem may be solved repeatedly, increasing resolution 
until the result converges, may also be employed here. 

Referring to FIG. 2 as an illustrative and non-limiting 
example, a grid size of 0.25 m may be deemed appropriate. 
This grid resolution is achieved by partitioning the faces of 
a 2 m cube as shown (into 64 identically-sized smaller 
cubes), or by modeling 0.25 m cubes with only 1 grid cell 
per face. The optimal choice will, of course, be problem 
specific, and involves trade offs between computational 
efficiency, error tolerance, and data requirements as dis- 
cussed in detail below. In any event, the size of cube chosen 
should be small enough to allow the model to approximate 
the physical shapes found in the modeled environment. 

Similar considerations apply to the determination of the 
directional resolution of the model. Monte Carlo simulations 
of radiative transfer generally partition a directional sphere 
onto discreet regions referred to as averaging quads or 
directional bins, reflecting the fact that the simulations do 
not keep track of the direction of each photon to pass the 
detection point, but instead simply increment the count of 
the appropriate directional bin. Customarily, the directional 
sphere is partitioned into quads bounded by circular arcs of 
constant polar angles (latitude) and circular arcs of constant 
azimuthal angles (longitude). The longitude lines are gen- 
erally evenly spaced, dividing the sphere into equally sized 
lines. The latitude lines may be drawn so that for all latitude 
lines, the difference in polar angles between any two adja- 
cent latitude lines is constant, or so that the difference in the 
cosine of the polar angles between any two adjacent latitude 
lines is constant. In either case, the triangular regions at the 
poles are often combined to form polar caps. 

According to the present invention, the former approach 
is preferred such that the triangular sections of the sphere are 
not combined into polar caps. This procedure results in bins 
in which the subtended solid angle varies depending on the 
polar angle of the bin. Bins near the sphere equator are 
larger, and this may be expected to result in collection of 
more “counts” than bins at the polar regions, implying that 
for any given Monte Carlo run, the values averaged over the 
equatorial bins will be more reliable than those averaged 
bins closer to the poles. This situation may be remedied by 
employing the latter method to draw latitude lines, but then 
the polar quads may be elongated compared to those nearer 
the Equator. Thus, two numbers are required to identify a 
directional bin; one to identify the polar zone, the theta bin 
number, and another to identity the azimuthal zone angle, or 
phi bin number. Thus, the raw output value for a face f, grid 
cell c, theta bin t, and phi bin p, denoted as O r (fc,t,p) is the 
total number of photons passing through face f, jgrid cell c, 
and the directional bin formed by the intersection of theta 
bin t and phi bin p. 

Once the spatial and directional resolutions are decided 
upon the memory requirements of the model are preferably 
checked to confirm that they remain tractable, though this is 
not required according to one aspect of the present inven- 
tion. One example is the memory requirements for the 
Monte Carlo modeling of the RF for the illustrative cube of 
FIG. 2. One copy of the RF will have ((6xl6x36x 
72)/2) 2 =l 5,479,341,056 data values, each of which requir- 
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ing 8 bytes for double precision calculations, which may 
exceed the performance abilities of many computers. 

Because the hybrid model approach can be tested on 
1 -dimensional and 2 -dimensional environments without loss 
of generality, the size of the RF’s can be reduced substan- 5 
tially. For example, if a 1 -dimensional environment is to be 
modeled at the resolution shown in FIG. 2, one method of 
constructing the model is to stack 3 -dimensional cubes in a 
column as shown in FIG. 4. Then, the output from the front 
face of the cube may serve as the input to the back face, and to 
vice versa. The output from the left face of the cube serves 
as the input to the right face, and vice versa. This is referred 
to as wrapping. Similarly, the output from the top face of a 
cube serves as the input to the bottom face of the cube above 
it, and vice versa. Even if all the cubes represent optically 15 
identical parcels, the full Monte Carlo model of the RF will 
be required. Then, for the iterative phase of the model, one 
copy of the RF and one copy of the output state for each cube 
in the column are required. This method of modeling 1 -di- 
mensional scenarios requires a significant amount of data 20 
storage and computation. 

According to the present invention, some or all of the 
wrapping may be done in the Monte Carlo modeling of the 
cube. For instance, assume the coordinate system local to the 
1 -dimensional environment to be modeled is situated such 25 
that the optical parameters vary in the z direction and are 
constant in the x and y directions. Wrapping the Monte Carlo 
model in the y direction is equivalent to modeling, instead 
of a cube, an infinitely long bar with square cross section as 
shown in FIG. 5. The number of faces is therefore reduced 30 
from 6 to 4. And because the square grid cells become long 
strips, their number is reduced from 16 to 4. The square 
plates used for boundaries also become strips. If the direc- 
tional resolution remains the same as before, the size of the 
RF for the bar is ((4x4x36 72)/2) 2 =429,981 ,696 data values, 35 
which is less than 3 percent of the former size discussed 
above. 

These bars may be stacked vertically and the outputs in 
the x direction wrapped to model the 1 -dimensional envi- 
ronment. To model a 2-dimensional scenario, an array of 40 
bars may be formed. A 3 -dimensional model will require the 
original 3 -dimensional cube to be modeled. Additionally, the 
Monte Carlo modeling of the bars is also performed in 
3 -dimensions. 

45 

ILLUSTRATIVE EXAMPLES AND TESTING OF 
THE EXAMPLES 

Two model environments may are next presented to 
illustrate the methods and advantages of the present inven- 50 
tion. A first example is a simple one-dimensional model 
representing a uniform water column, depicted in FIG. 6. 
The second example is a two-dimensional scenario in which 
a long barge floats over a uniform water column, as shown 
in FIG. 7. In the second example the bottom and sides of the 55 
barge were modeled as Lamberian reflectors with a reflec- 
tance of 0.1, and the right and left sides of the two- 
dimensional environment were wrapped in the above-de- 
scribed manner, creating a model environment representing 
a series of long parallel barges, each 10 m wide and spaced 60 
20 meters apart. In both examples, the bottom and air/sea 
interface are flat and the water is 1 0 m deep. The sky for both 
cases consists of a direct beam at a zenith angle of 30° 
embedded in a uniform, difluse background sky. The diffuse 
sky provides 21% of the downward irradiance. 65 

Both environments were constructed using the 2 -dimen- 
sional bars described above, each having 8 grid strips per 
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face. The directional sphere was divided into 1 8 theta bins 
and 20 phi bins. The size of the RF for these bars is 
((4x8x1 8x20) 2 =33, 177,600 data values, which, for double 
precision accuracy, requires about 265 MByte of data stor- 
age. 

For each environment, two water types were modeled: 
one representing clear, offshore water, the other more turbid 
water characteristic of ports and estuaries, such as Tampa 
Bay. The optical parameters were set to match typical 
chlorophyll concentrations of 0.3 mg/m 3 and 5.0 mg/m 3 , 
respectively. These readings are tabulated in Table 1 , below. 
Additionally, both cases were simulated at a wavelength of 
532 nm and the bottom reflectance was set to 0.15 for the 
clearer water and 0.05 for the more turbid estuarine water. 


TABLE 1 


Optical Parameters For Two Water Types, 

Wavelength = 532 nm 

Chi - 0.3 mg/m 3 Chi - 5.0 mg/m 3 

Absorption (water) 

0.05172 nr 1 

0.05172 nr 1 

Scattering (water) 

0.00218 nr 1 

0.00218 nr 1 

Absorption (particles and 
Colored Dissolved 

0.01304 m" 1 

0.1150 m -1 

Organic Matter) 
Scattering (particles) 

0.1448 nr 1 

0.8391 m -1 

Bottom reflectance 

0.15 

0.05 

Barge reflectance 

0.10 

0.10 


One-Dimensional Example 

Downward irradiance calculated by the hybrid model 
compares well with the Hydrolight results as shown in FIG. 
8. The RMS difference between the two models for the clear 
water is 0.98%. For the turbid water, the difference is 2.36%. 
In both cases, the downward irradiance transmittance of the 
hybrid model is greater than that of Hydrolight. 

Upward irradiance calculated by the hybrid model is 
comparable, but the difference in the two models is slightly 
increased: 1.28% RMS for the clear water, 3.63% for the 
turbid water as depicted in FIG. 9. Part of the difference in 
upward irradiance can be attributed to the difference in 
downward irradiance at the bottom. Even so, in this case the 
upward irradiance transmittance of the hybrid model is less 
than that of Hydrolight. Results shown in FIGS. 8 and 9 
imply that the hybrid model may transport radiance prefer- 
entially depending on the direction of transport and the 
orientation of the model element. In this case, the direct solar 
beam may be transported more effectively than the diffuse 
radiance reflected from the bottom. 

FIG. 10 shows the average cosines vs. the depth calcu- 
lated by the two models. The upward average cosine is the 
average value of the polar angle of all the photons contrib- 
uting to the upward irradiance at the given depth. The 
downward average cosine is defined analogously. The dif- 
ferences between the two models range from 1 .40% RMS to 
2.48% RMS, indicating that the subsurface radiance distri- 
butions generated by the two models are practically equiva- 
lent. Differences in the downward average cosines just 
beneath the surface may be related to differences in the way 
the two methods account for transmission and reflection at 
the interface. 

The total upward radiance just above the surface in the 
half plane perpendicular to the principal plane is shown for 
both models and both water types in FIG. 11 . The difference 
between the two models for the clear water is 7.36% RMS. 
For the turbid water, the difference is 8.94% RMS. 
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Several effects may contribute to the disparity between Two-Dimensional Example 

the two radiance distributions. First, in the present example 

the more than 33 million values defining the RF for the bar FIG. 12 shows the normalized downward irradiance 
elements were derived in the Monte Carlo simulation from beneath and adjacent to a barge afloat in the clear water type, 

no more than 100 million weighted photon packages. The 5 it i s important to note that vertical and horizontal axes are 
calculation of irradiances and average cosines requires the plotted to different scales. Data fortheseplots (FIGS. 12-15) 

combination of the output values for many counting bins, was calculated at an eighth-meter resolution, then smoothed 

thus smoothing occurs of the rough data resulting from the t0 a spatial resolution of one meter for plotting. The vertical 

dearth of sample trajectories. Calculation of the radiance ^ lines just below the bottom comers of the barge are plotting 
distributions combines fewer output values; so increased artifacts caused by the combination of this low spatial 
variance is expected. For example, the total upward irradi- resolution and the high horizontal gradient of downward 
ance just above the surface differs by only 2 . 25 % for the two irradiance at the edges of the barge. The shadow that extends 

models for the clear water case, and this difference would above the sea surface on the left side of the barge is clearly 
likely be reduced if more photon packages were used to 15 shown, as is the sharp horizontal gradient in downward 
determine the RF for the bar elements. But even if this value irradiance at the lower right comer of the barge. The 
is considered to be a systematic bias between the two contours indicate that a submersible vehicle moving beneath 
models, simply increasing the number of photon packages t h e barge at a depth of four meters would encounter changes 
used to determine the RF reduces the random errors in the j n downward irradiance of nearly three orders of magnitude, 
radiance distributions to a comparable level. 20 This downward irradiance diffuses beneath the barge to fill 

Secondly, the differences shown in FIG. 11 represent the in the shadow, but at the same time, the downward irradiance 

accumulation of all the errors in the hybrid model because at the bottom in the deepest part of the barge shadow is more 

the upward radiance distribution above surface relies on the than an order of magnitude less than in the open spaces 

transmission of the downward light from the source, through adjacent to the barge. Curvature of the downward irradiance 

the air-sea interface, to the bottom of the water column, and 25 contours extends to the edges of the plot, indicating that the 
then back. Therefore, the differences in the radiance distri- downward irradiance in the open areas between the barges 

butions effectively represent the upper limit of the model is perturbed by the presence of the baiges and shadows, 

error. Even so, the hybrid model results mimic the shape of FIG. 13 shows the upward irradiance modeled for the 
the Hydrolight-produced radiance distribution quite well. ^ same conditions prevailing for FIG. 12. Because most of the 
Finally, there are fundamental differences in the structure upward irradiance originates from reflection off the Lam- 
of the two models that make it difficult to formulate perfectly bertian bottom and is more diffuse than the downward 
equivalent problems. Hydrolight uses the polar cap irradiance, the structure in the upward irradiance seen at the 
described above to combine all the bins above a zenith angle bottom dissipates with height. Unlike the downward irradi- 
of 5 °, and then partitions the remainder of the hemisphere 35 ance field, the upward irradiance is more uniform just 
down to 85 ° into 10 ° zonal bins centered at zenith angles of beneath the barge than it is just above the bottom. It is noted 
10 °, 20 °, etc. The hybrid model, at the directional resolution that the upward irradiance curvature contours to the right of 
used for this work, partitions the hemisphere into 10 ° zonal the barge. The upward irradiance reaches a maximum for 
bins centered at zenith angles of 5 °, 10 °, etc. Thus, Hydro- each depth at about 28-29 meters on the horizontal scale 
light spreads the direct solar beam across a bin spanning 4Q indicating that for this combination of model geometry and 
from 25 ° to 35°, whereas the hybrid model spreads the beam optical parameters, the series of barges are spaced so closely 
across a bin spanning from 20 ° to 30 °. This may account for that the upward irradiance at any point more than a few 
the higher downward average cosine at the surface generated meters from the bottom is affected by the shadows of at least 
by the hybrid model, and also for the slightly higher down- two barges. 

ward irradiance transmittance. 45 FIGS. 14 and 15 show the results of the hybrid method 

There are also differences in the way the two models when the turbid water parameters are modeled in the same 
handle sub- surface internal reflection. The largest disparities geometry as for FIGS. 12 and 13. In this case, the downward 

in the radiance distributions shown in FIG. 11 occur for irradiance reflected from the bottom of the barge is about 

zenith angles greater than 75 °, the portion of the distribution four orders of magnitude less than the total downward 

most affected by the subsurface internal -reflection algo- 50 irradiance at the surface. The shape of the downward irra- 
rithm. When data points corresponding to zenith angles diance field (FIG. 14) is similar to that shown in FIG. 12, 

greater man 75 ° are not considered, the RMS differences except that the contours of FIG. 14 are closer together 

between radiance distributions generated by the two models because of the increased attenuation of the turbid water. At 
for the clear and turbid water types drops from 7 . 36 % RMS the edges of the plot, the downward irradiance contours are 
to 5 . 22 % RMS and from 8 . 94 % to 7 . 35 %, respectively. 55 approximately linear and parallel to the bottom, an indica- 
Thus, so even within the Hydrolight model, use of finite bin tion that the downward irradiance in this region is not 
sizes that don’t precisely match and center about the geom- affected by the presence of the barges, 

etry of the source and viewing angles will cause small The upward irradiance field (FIG. 15), however, is quite 

discrepancies. different than that shown in FIG. 13. The vertical contours 

The above-surface radiance distributions should be axi- 60 emanating from the bottom of the barge in FIG. 15 show that 
ally symmetric when the sky conditions are changed for the the gradient of upward irradiance just beneath the barge is 
model scenario so that the sun is at zenith with a uniform horizontal. The bottom of the barge is illuminated solely 
background sky. When the Hydrolight model is executed from infilling of irradiance from the open areas between the 
with a zenith sun, the resulting above-surface radiance barges, not by reflection from the bottom. Unlike the clear 
distributions in the half-planes at 45 ° and 90 ° to the principle 65 water case, the upward irradiance distribution in the center 
plane differ by 1 . 28 % RMS and 1 . 19 % RMS for the clear off the region between the barges appears to be unperturbed 
and turbid waters, respectively. by the barges or their shadows. For both water types, the 
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effect of the reflection of downward irradiance from the 
bright side of the barge is shown clearly by the upward 
irradiance contours. 

Applications of the Methods of the Present 5 

Invention 

Optical sensors typically work within limited ranges of 
radiance. In order to optimize sensor performance, it is 
necessary to define the range of radiance within which the 10 
sensor is expected to function. Special operating techniques 
for sensors which are deployed in widely varying light fields 
may be necessary. For example, the light field beneath the 
barge as shown in FIG. 14 is more than three orders of 
magnitude lower than it is at the sea surface. To compensate 15 
for such a range, an automatic gain control (AGC) may be 
used to allow a single sensor or camera on an autonomous 
underwater vehicle (AUV) or remotely operated vehicle 
(ROV) to be effective in both environments. Furthermore, if 
the field of view of a video camera included both the right 20 
side and bottom of the barge, as shown in FIG. 12, the 
camera will view portions of the barge that contrast by more 
than two orders of magnitude. The automatic gain control 
for the sensor will then adjust to a value that is far too 
insensitive for use in discriminating objects on the bottom of 25 
the barge. Only when the side of the barge is excluded from 
the field of view will the inspection of the barge bottom be 
feasible. Operational practices can be devised to avoid such 
problems, one method being the use of zoom lens adjust- 
ments to avoid bright pixels. 30 

Because of the more reflective bottom and less turbid 
waters typical of an offshore environment, almost three 
times as much natural light is available to illuminate the 
bottom of a ship for imaging by video camera or other 
conventional passive imaging sensors in offshore waters as 35 
is available in an inshore environment. On the other hand, if 
a laser line scanner is employed to inspect the bottom of a 
ship during daylight, higher contrast for observing the laser 
line would be available beneath a ship in inshore water. The 
hybrid model can produce the estimates of the ambient light 40 
field needed to determine the optimal location and instru- 
ment to be used for hull inspection. 

Inspection of ship hulls is time-consuming, and therefore 
expensive; the hull of a ship transports liquefied natural gas 
may span more than five acres. Even so, because the 45 
consequences of such a vessel exploding in port would be 
catastrophic, divers now take up to two weeks to inspect the 
hulls of such ships before allowing entry to ports such as 
Boston Harbor. Delays in making port due to the inspection 
process could be greatly reduced using optical systems 50 
presently being developed when mounted on AUVs and/or 
ROVs. The savings accrued by employing such systems to 
inspect a large vessel would be on the order of $100,000 to 
$250,000 for each day by which the inspection process was 
shortened. 55 

To prepare the operators of such vehicles for inspections 
beneath ships in various environments, a preview of the 
environment should be available in order to optimize the 
range and observational parameters for best sensor perfor- 
mance. The hybrid model can be run a priori to provide good 60 
estimates of the optical environments to be expected for a 
wide range of conditions. Then, given contemporaneous 
values of the optical qualities of the water and natural light 
conditions acquired by remote sensing or instrumented 
buoys, be vehicle operator may quickly extract the optimal 65 
sensor parameters from a look-up table. The parameters 
measured include water conductivity, temperature, sun 
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angle, cloudiness, depth, water turbidity and any other 
parameters, such as absorptional and scattering coefficients, 
as known to one of ordinary skill in the art that affect the 
transmission of light in the test medium. 

The instant invention has been described in terms of an 
underwater environment but it is within the scope of the 
invention that any other environment may also be modeled. 
This includes, but is not limited to, human, other mammalian 
or other types of body, including marine life, reptiles, etc. In 
addition, other environments, such as subaerial mountainous 
settings as known to those of skill in the art may also be 
modeled as long as they have methods of optical measure- 
ment applied to them. 

In summary, by use of the Monte Carlo and iterative 
techniques simulated solutions to the Radiative Transfer 
Equations in one-, two- and three-dimensions can be made. 
The optical environment to be modeled is first partitioned 
into contiguous elements. In the illustrative example 
described herein, the elements were cubes or infinitely long 
bars of square cross-section. For each type of element, the 
type being defined by the geometric shape of the element 
and the distribution of optically active elements within it, 
Monte Carlo techniques are used to determine the element’s 
optical response function or RF. Once the RF is calculated 
for every type of element in the region, the elements are 
assembled (conceptually), and the boundary conditions are 
imposed. Boundary conditions are implemented using spe- 
cialized elements with specified, constant output states to 
represent sources, and elements with analytically deter- 
mined response functions to represent reflective surfaces and 
Fresnel interfaces. The outputs of the source elements are 
diffused iteratively throughout the model region until the 
output states of all elements change negligibly with further 
iteration. One-dimensional results of the hybrid model com- 
pare well with the widely used Hydrolight program, with the 
largest differences seen in the above-water total upward 
radiance distributions. 

The advantages of the hybrid model are exemplified by 
simulation of the two-dimensional light field adjacent to and 
beneath a series of evenly spaced long barges. The model 
region was partitioned into long bars with one-meter square 
cross-sections. The sides of each bar were further divided 
into eight strips. Because the water in the model region was 
defined to be homogeneous, only two types of bar elements 
were necessary; one type each for the clear and turbid water 
simulations. The iterative part of the hybrid model method 
combined the bars and solved simultaneously for the radi- 
ance distribution at each strip on each face of every bar in 
the model region. In other words, the complete radiance 
distribution, with a directional resolution of 1 8 theta values 
for each of 20 phi values, was determined at about 5000 
locations throughout the model region using only one Monte 
Carlo simulation for each of the two water types comprising 
the model environment. 

By use of this hybrid model radiance estimates are 
sufficient to allow for sensor calibration as well as predict- 
ability of performance, as well as whether or not data 
collection is viable for any set of parameters. This enables 
ship hull and port facility inspection to be done by ROV and 
AUV systems, that otherwise would require more costly and 
time-consuming techniques. In addition, the methods of the 
present invention apply to systems wherein the instrumen- 
tation is located on a spacecraft or on any type of aircraft 
system, including piloted, autonomous, or semi -autonomous 
systems. Combinations of these systems are also contem- 
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plated within the scope of the instant invention, so that 
simultaneous monitoring from a plurality of locations is also 
possible. 

Modification and variation can be made to the disclosed 
embodiment of the instant invention without departing from 5 
the scope of the invention as described. Those skilled in the 
art will appreciate that the applications of the present 
invention herein are varied, and that the invention is 
described in the preferred embodiment. Accordingly, addi- 
tions and modifications can be made without departing from to 
the principles of the invention. Particularly with respect to 
the claims it should be understood that changes may be made 
without departing from the essence of this invention. In this 
regard it is intended that such changes would still fall within 
the scope of the present invention. Therefore, this invention 15 
is not limited to the particular embodiments disclosed, but is 
intended to cover modifications within the spirit and scope 
of the present invention as defined in the appended claims. 

What is claimed is: 

1. A method of determining a radiance field in an optical 20 
environment, said method comprising the steps of: 

generating a model of the optical environment; 

partitioning the modeled, optical environment into dis- 
creet, contiguous regions; and 

applying Monte Carlo techniques to determine an optical 25 
response function for each region for which the optical 
response function cannot be analytically defined; and 

reporting the optical response function determined by 
application of the Monte Carlo techniques. 

2. A method as set forth in claim 1 further comprising the 30 
step of combining the optical response function for each 
region. 

3. A method as set forth in claim 2 wherein the step of 
combining the optical response function for each region 
comprises utilizing an output of one region as an input to a 35 
contiguous region. 

4. A method as set forth in claim 2 further comprising the 

step of applying iterative relaxation techniques to determine 
the radiance field on a boundary and throughout the mod- 
eled, optical environment. 40 

5. A method as set forth in claim 1 wherein the step of 
applying Monte Carlo techniques to each region comprises 
solving the Radiative Transfer Equation (RTE) for each 
region to obtain an end solution. 

6. A method as set forth in claim 5 further comprising the 45 
step of calibrating data from an optical sensor in response to 
the end solution. 

7. A method as set forth in claim 6 wherein the step of 
calibrating data comprises adjusting gain of the optical 
sensor in response to the end solution. 
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8. A method as set forth in claim 6 wherein the optical 
sensor is located on at least one of a remotely operated 
vehicle (ROV), an autonomous underwater vehicle (AUV), 
a piloted aircraft, an autonomous aircraft, a semi-autono- 
mous aircraft, a piloted spacecraft, an autonomous space- 
craft, and a semi -autonomous spacecraft. 

9. A method as set forth in claim 5 further comprising the 
step of validating data from an optical sensor in response to 
the end solution. 

10. A method as set forth in claim 9 wherein the optical 
sensor is located on at least one of a remotely operated 
vehicle (ROV), an autonomous underwater vehicle (AUV), 
a piloted aircraft, an autonomous aircraft, a semi-autono- 
mous aircraft, a piloted spacecraft, an autonomous space- 
craft, and a semi -autonomous spacecraft. 

11 . A method as set forth in claim 5 further comprising the 
step of determining an optical sensor to deploy in response 
to the end solution. 

12. A method as set forth in claim 5 further comprising the 
step of adjusting an optical sensor in response to the end 
solution prior to deployment of the optical sensor. 

13. A method as set forth in claim 1 wherein the optical 
environment is underwater. 

14. A method as set forth in claim 1 wherein the optical 
environment is subaerial. 

15. A program product for determining a radiance field in 
an optical environment wherein program product code is 
stored on computer readable media and comprises: 

a modeling code stored on the computer readable media 
for generating a model of the optical environment; 

a partitioning code stored on the computer readable media 
for partitioning the modeled, optical environment into 
discreet, contiguous regions; and 

an application code stored on the computer readable 
media for applying Monte Carlo techniques to deter- 
mine an optical response function for each region for 
which the optical response function cannot be analyti- 
cally defined. 

16. A program product as set forth in claim 15 further 
comprising a combination code stored on the computer 
readable media for combining the optical response function 
for each region. 

17. A program product as set forth in claim 16 further 
comprising a relaxation code stored on the computer read- 
able media for applying iterative relaxation techniques to 
determine the radiance field on a boundary and throughout 
the modeled, optical environment. 





